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Within a cosmological context, we study the behaviour of collisionless particles in the weak field 
approximation to General Relativity, allowing for large gradients of the fields and relativistic veloci¬ 
ties for the particles. We consider a spherically symmetric setup such that high resolution simulations 
are possible with minimal computational resources. We test our formalism by comparing it to two 
exact solutions: the Schwarzschild solution and the Lemaitre-Tolman-Bondi model. In order to 
make the comparison we consider redshifts and lensing angles of photons passing through the sim¬ 
ulation. These are both observable quantities and hence are gauge independent. We demonstrate 
that our scheme is more accurate than a Newtonian scheme, correctly reproducing the leading-order 
post-Newtonian correction. In addition, our setup is able to handle shell-crossings, which is not 
possible within a fluid model. Furthermore, by introducing angular momentum, we find configu¬ 
rations corresponding to bound objects which may prove useful for numerical studies of the effects 
of modified gravity, dynamical dark energy models or even compact bound objects within General 
Relativity. 


I. INTRODUCTION 

Recent results from the Planck mission [1], BOSS [2-4], 
WiggleZ survey [5], CFHTlenS [6] and SNLS [7] have con¬ 
solidated the ACDM concordance model of cosmology as 
providing a very good fit to observations. However, this 
model is characterized by two semi-phenomenological in¬ 
gredients - cold dark matter (CDM) and a cosmological 
constant (A) - whose true nature still needs to be de¬ 
termined at the fundamental level. With the reach of 
linear analysis now being nearly exhausted, phenomena 
at nonlinear scales can help to make progress. On the 
observational side, large surveys such as Euclid [8], DES 
[9], LSST [10] and SKA [11] will make significant progress 
analysing these non-linear scales. This puts the onus on 
the theoretical side to understand precisely what we ex¬ 
pect these surveys to see. Due to the non-linearity, nu¬ 
merical simulations will be a necessary tool to probe this 
regime. 

The N-body codes used for the study of cosmic large 
scale structure normally employ Newton’s law of gravi¬ 
tation. One expects that this approximation works well 
as long as perturbations are generated by nonrelativis- 
tic matter only. This is true if dark energy is indeed 
a cosmological constant and dark matter is some heavy 
fundamental particle (like in the WIMP scenario). How¬ 
ever, since these facts are not established it seems that 
by using the Newtonian approximation we are unable to 
access a viable part of model space. 

In fact, even within the realm of known physics this ap- 
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proximation will break down due to the existence of very 
light, but still massive, neutrinos. In principle, the ini¬ 
tial conditions of simulations can be set late enough that 
neutrinos have already become non-relativistic; however, 
this can be so late that the cold dark matter has already 
begun to cluster significantly. Therefore, relativistic ef¬ 
fects are already important in order to rigorously model 
the effects of neutrino masses in cosmology. 

In an effort to address these shortcomings, a relativis¬ 
tic framework for N-body simulations has recently been 
developed [12, 13]. This framework is based on a weak- 
field expansion of Einstein’s equations, similar to the one 
proposed in [14, 15]. It does not require a particular 
form of stress-energy and relies solely on the assumption 
that gravitational fields are weak, at least at large scales. 
Therefore, it is applicable to a much larger set of mod¬ 
els, including hot dark matter [16, 17] and many types of 
dynamical dark energy [18]. 

Before investing significant computational resources in 
order to do a full-scale cosmological simulation it is in¬ 
teresting to study the relativistic effects in a simplified 
setup. Here we will consider the case of a single, iso¬ 
lated, spherically symmetric structure which could, for 
instance, be a model for a cosmological void or a galaxy 
cluster. The idealization to exact spherical symmetry 
drastically reduces computational requirements, allowing 
high-resolution simulations to be carried out at negligible 
cost. Furthermore, the numerical scheme can be thor¬ 
oughly verified by comparing to several known exact so¬ 
lutions. When comparing to exact solutions, structures 
can also be allowed to evolve into regimes where metric 
perturbations do become large and the framework breaks 
down, allowing us to probe the boundaries of where the 
framework can and cannot be trusted. 

Our approach is in some sense complementary to ex¬ 
isting methods for the numerical solution of Einstein’s 
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equations. For instance, the BSSN formalism [19-21] can 
probe the strong field regime, but existing implementa¬ 
tions rely on a fluid description for matter. Our N-body 
method, on the other hand, allows us to study matter 
configurations with highly nontrivial phase space distri¬ 
butions. 

In section II we introduce the relativistic framework 
and study some simple spherically symmetric setups. 
We first consider a Schwarzschild solution to confirm 
that the relativistic potentials are calculated accurately 
in vacuum. We then add nonrelativistic matter and 
compare our simulations to the exact Lemaitre-Tolman- 
Bondi models which describe spherically symmetric so¬ 
lutions with a dust fluid. In order to avoid gauge issues, 
we construct several physical observables which can be 
compared without ambiguity. We note that the fluid 
solutions break down at the formation of caustics, but 
our relativistic framework remains valid and can thus 
probe settings beyond the fluid approximation. Without 
support from pressure or angular motion, overdensities 
tend to collapse quickly and can not easily form stable 
bound objects. In section III, we propose a way to in¬ 
troduce angular motion without breaking spherical sym¬ 
metry. This is achieved by arranging the motion of the 
particles such that they all individually have angular mo¬ 
mentum, but the total angular momentum of the system 
remains zero. We demonstrate that one can And config¬ 
urations corresponding to bound objects. Such configu¬ 
rations may be useful laboratories to study the effects of 
modified gravity, dynamical dark energy models, or even 
the early stages of the formation of primordial blackholes 
[22, 23], or ultra compact mini-haloes [24-26], within or¬ 
dinary gravity. 

II. THE MODEL 

The perturbed Friedmann-Lemaitre-Robertson- 
Walker (FLRW) metric, in spherical coordinates and 
longitudinal gauge, is: 

ds^ = —a^(T) [1 -I- 2'I'(r, r)] dr^ 

-I- a^(r) [1 — 2<i)(r, r)] [dr^ -|- r^dfl^] , (1) 

where r is the conformal time, a is the scale factor and 
we impose spherical symmetry of the perturbations by 
requiring that the Bardeen potentials <i)(r, r) and r) 
depend only on the radial coordinate and time. We 
have also assumed a spatially flat background although it 
would be easy to generalize our model to allow for open 
or closed geometries. 

We examine this metric in the regime where gravita¬ 
tional fields are weak. In other words, we are interested 
in perturbations caused by structures that remain much 
larger than their Schwarzschild radius. Such a weak-fleld 
setting allows for a systematic expansion of the various 
equations of motion (including Einstein’s held equations) 
in terms of metric perturbation variables. We will follow 


an approach studied in [12] which takes into account the 
most important relativistic terms. 

This approach can be summarized as follows: first, all 
equations are expanded in terms of the metric pertur¬ 
bations - in our case <i> and T - and all terms up to 
first order are kept without distinction. At higher orders, 
however, one only wants to keep the most relevant terms. 
Noting that linear perturbation theory is accurate on the 
largest scales (close to or beyond the horizon) the only 
higher order terms that we will keep are those which may 
become large at small scales. These terms will be those 
with two spatial derivatives,^ since a derivative will ef¬ 
fectively multiply a term by an inverse power of a length 
scale. To arrive at a tractable set of equations that still 
contains the most important relativistic corrections we 
will therefore add all second order terms with two spatial 
derivatives and no terms of any higher order. Although 
there are scenarios where terms of higher than quadratic 
order can dominate over the linear terms (e.g. <i> y if 
^,ij > ^ij /‘I’) these higher order terms will always be sub¬ 
dominant to the largest quadratic order terms that we do 
include. Further details on this approximation scheme 
can be found in [12]. 

It is important to emphasise that any perturbations of 
the stress-energy tensor, including momenta, are allowed 
to be arbitrarily large. The perturbative expansion is 
only carried out in terms of gravitational fields and we 
make no assumptions about other perturbations. For 
instance, our solar system perfectly fits into this scheme 
since the gravitational held of the sun remains well within 
the weak-fleld regime, despite the fact that its density is 
some thirty orders of magnitude larger than the mean 
density of the Universe. 

Using the “time-time” component of Einstein’s equa¬ 
tions, Gq = SttGTq , we obtain an equation for the metric 
perturbations: 

$ + - 3H2($ -x)+ ^(^.r)" 

r 2 

= -AnGa'^il - 44>)JT° , (2) 

where commas denote partial derivatives with respect to 
r or r. Note that, in a spherical coordinate system as 
the one used here, second spatial derivatives can give 
rise to terms like 4* ^/r. We will treat these terms like 
second derivatives in our expansion scheme, i.e. a fac¬ 
tor 1/r will effectively be counted like a spatial deriva¬ 
tive. We also introduced the conformal Hubble parame¬ 
ter 7i = dlna/dr and the difference of the potentials as 
X = 4) — 'k. On the right-hand side, STq stands for the 
perturbations of the stress-energy tensor, STq =Tq—Tq. 
We will only consider contributions from massive parti¬ 
cles (e.g. cold dark matter). The background model is 


^ Note that Einstein’s equations are second order differential equa- 
tions, therefore no terms will have more than two derivatives. 
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governed by the Friedmann equation 

n^ = ~a^fS. (3) 

Another equation comes from the traceless part of the 
“space-space” components of Einstein’s equations, G* — 
{'^j ~ reads: 

X,rr — ^X,r + X^r + + 2 ~ X) 

= 127rGa2(l - 2x)n„, (4) 

where 11^.^ is the radial component of the anisotropic 
stress, defined for a general coordinate system as: 

n., = . (5) 

Latin indices denote spatial coordinates only. As a con¬ 
sequence of spherical symmetry the anisotropic stress is 
purely longitudinal in our setting. 

The stress-energy tensor is derived by varying the ac¬ 
tion of an ensemble of massive point particles with re¬ 
spect to 5g^u (see e.g. equation (2) of [13]). This gives: 




(5(3)(x-x(„)) 


^ I 9a(3 


V^9 
d^fn) ’ 

dr dr 


- 1/2 


dx'f„. da; 


(") 


dr dr 


( 6 ) 


where we sum the contributions of n particles with 
masses m(„) and spatial positions and Greek in¬ 

dices run over all four coordinates of space-time. 

For a particle moving in the radial direction we can 
define a momentum 


P = 




1-^24-- (1-2$) (f^) 


( 7 ) 


which is the proper relativistic momentum as measured 
in a Gaussian orthonormal coordinate frame aligned with 
our foliation of spacetime^. The motivation for this is 
that it allows us to derive an expression for the stress- 
energy tensor that is valid (within the bounds of our ap¬ 
proximation scheme) even when the particles have arbi¬ 
trarily high velocities. In particular. 


^ Explicitly, a Gaussian orthonormal coordinate frame is given by 
a set of orthonormal basis vectors Cq , e'^, = — 1, 

= 0, = 5ij with orthogonal to the spacelike 

hypersurface. The metric in the coordinates defined by this basis 
locally looks like the Minkowski metric. The momentum p de¬ 
fined in eq. (7) is simply the spatial component of the covariant 
4-momentum in that coordinate system. If we align with the 
radial direction, i.e. oc <5^ , we can write p = where 

denotes the covariant 4-velocity. 


and 




2 1 + 3$ 

3 47rr^a^ 




r(n)) 



(9) 


As we restrict our solutions to spherical symmetry we 
can imagine collections of particles as representing spher¬ 
ically symmetric shells with only radial positions. This 
is because symmetry also requires that the particle dis¬ 
tribution function is independent of angular position. 
We therefore simply dropped the angular coordinates 
from above expressions and defined the masses such that 
m(„) / (47rr^^^) is the surface mass density (in coordinate 
space) of the shell with label n. Thus, each shell ac¬ 
counts for all particles at given radius with given 
radial momentum and we need only sum over shells. 
Note that, for descriptive ease, from here onwards we 
will refer to the individual shells as the “particles” of our 
simulations. 

In order to evolve the particle positions one can invert 
eq. (7), 


dr 

dr 


P 

\J m? 


(l-h$-k$). 


The geodesic equation for massive particles, 
ds2 ^ ’'P ds ds ’ 


( 10 ) 


( 11 ) 


finally determines the evolution of the momenta as 


dp 

dr 


— {H — $,t)p — -I- p2 _ 


( 12 ) 


Our aim is to study numerical solutions to the above 
system of equations. To this end, we adopt a particle- 
mesh (PM) scheme as used in many cosmological N-body 
simulations. The “mesh” part of the scheme takes care 
of the evolution of fields such as $ or y. All fields are 
represented approximately by sampling their values on a 
discrete set of points, hereafter referred to as the “grid”. 
The field equations (2), (4) are solved on the grid by re¬ 
placing the differential operators by finite-difference ver¬ 
sions thereof. 

The “particle” part of the PM scheme, on the other 
hand, takes care of the evolution of the particle ensemble. 
The phase-space of fundamental particles is sampled by 
a much smaller number of N-body particles which can 
be viewed as discrete elements of phase-space. Hereafter, 
the term “particle” usually refers to the latter notion. 
Positions and momenta of particles are real-valued (i.e. 
they can exist in arbitrary positions between grid points) 
and the geodesic equation is solved by interpolating field- 
dependent quantities such as to the particle positions. 

Vice versa, a so-called particle-to-mesh projection is re¬ 
quired to construct the stress-energy tensor (whose com¬ 
ponents are treated like a field) from the particle ensem¬ 
ble. This is achieved by replacing J(r—r(„)) —> w(r—r(„)) 
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in eqs. (8), (9), where w is a weight function which de¬ 
pends on the projection method. We use the so-called 
“triangular-shaped particle” (TSP) method where w is 
constructed using a piecewise linear (triangle-shaped, 
hence the name) function of the separation. Some de¬ 
tails on the projection and interpolation methods can be 
found in appendix C. 

In the following subsections, in order to validate the 
numerical scheme, we will compare simulations to two 
well-known exact solutions of Einstein’s equations. 


A. The Schwarzschild Solution 


The Schwarzschild metric describes the spherically 
symmetric vacuum solution around a central mass con¬ 
centration. Within the context of our simulations this 
metric is suitable for regions void of particles, and can 
hence be used to test the implementation of the field 
equations independently of the particle evolution. 

In order to obtain explicit expressions for the two 
Bardeen potentials, it is useful to write the Schwarzschild 
solution in so-called “isotropic coordinates” [27], 


ds2 = _I^^- + (13) 

(1+^)2 V 4r; L J > i ^ 


where rg = 2GM denotes the Schwarzschild radius. 

As long as the mass M is distributed over a central re¬ 
gion much larger than rg, the exterior Schwarzschild so¬ 
lution can be viewed as a perturbation around Minkowski 
space. Within our simulations, such a background is de¬ 
scribed by Tq = 0 and hence H = 0. We can therefore set 
0=1 and T = t. In order to obtain a numerical solution 
we set up a homogeneous ball of particles (much larger 
than its Schwarzschild radius) in the center of an other¬ 
wise empty simulation volume. The Bardeen potentials 
'k and 4) outside of the ball are independent of time, as 
guaranteed by Birkhoff’s theorem. Their behavior in the 
weak-field regime is given by the large-r expansion of the 
above exact metric. For r ^ rg we have 

l + 2^{r) =1-!^ + ^ + ..., (14a) 

1 - 2$(r) =1+^ + Y4 + ... . (14b) 
r 8r^ 

In Figure 1 we show some results for 4> and x ob¬ 
tained with our numerical scheme and compare them to 
the corresponding analytic results obtained from the ex¬ 
act solution. We only consider the vacuum region out¬ 
side the central mass concentration. Evidently, our nu¬ 
merical scheme accurately accounts for the leading-order 
post-Newtonian corrections and is therefore one order (in 
post-Newtonian counting) better than a purely Newton¬ 
ian scheme. Using the results of our simulation it would 
be possible, for instance, to get an accurate prediction 
for the advance of the perihelion of Mercury. 
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Figure 1: (Color online) Top: numerical results for 4> (green) 
and X = <& — 4^ (red) as a function of r/rg. They are in ex¬ 
cellent agreement with the exact result, shown as dashed and 
dot-dashed lines, respectively. Bottom: relative error of the 
numerical values of $ with respect to the exact result (in 
green). For comparison, the relative error of the Newton¬ 
ian approximation, = —GM/r, is shown as dotted line (in 
black). 


B. The Lemaitre-Tolman-Bondi Solution 

If spacetime is filled with a dust fluid then one can 
construct a spherically symmetric class of exact paramet¬ 
ric solutions known as Lemaitre-Tolman-Bondi (LTB) 
models. For our simulation this is equivalent to requir¬ 
ing that particles occupying identical space-time points 
also have identical velocities, or more precisely, having 
a phase space distributon function which, at each space- 
time point, is a single Dirac delta-distribution in velocity 
space. Being exact solutions, these models do not require 
the density to be nearly homogeneous, allowing the study 
of strongly non-perturbative settings. However, since the 
construction makes use of a comoving synchronous co¬ 
ordinate system in an essential way, it is impossible to 
extend these solutions beyond the point where particle 
trajectories cross and the phase space distribution func¬ 
tion loses its simple delta-distribution character. In this 
coordinate system, the LTB line element reads 

+ f+2F(r) ^"' + 

In order to compare the LTB solutions to our numer¬ 
ical calculations we choose initial conditions such that 
the density perturbation is linear. In the perturbative 
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regime we can easily work out the gauge transformation 
between the synchronous comoving coordinates used to 
parametrize the LTB solution and the coordinates used 
in our framework, which are related to the longitudinal 
gauge. Details can be found in appendix A. Setting iden¬ 
tical initial conditions in this way we can make a com¬ 
parison based on identical physical situations. 

Figures 2 and 3 show some simulation results. We plot¬ 
ted the evolution of the density contrast, momentum of 
particles p, scalar perturbation and the difference of the 
two potentials y as a function of comoving radius r. The 
first set of figures portrays the collapse of a spherically 
symmetric overdensity and the second set shows the ex¬ 
pansion of a spherically symmetric void. The initial den¬ 
sities were set as “compensated tophat profiles”, where a 
central region of constant density contrast is surrounded 
by a second layer with constant density contrast of oppo¬ 
site sign such that the entire region can be matched onto 
a homogeneous FLRW exterior solution. In both cases 
we find X to be proportional to ^ and negative. 

In the expansion of a void, we can also observe a “shell¬ 
crossing” which happens when a set of particles moves 
outwards faster than particles that were initially at a 
larger radius. As mentioned before, this cannot be rep¬ 
resented with an LTB solution - it becomes singular as 
soon as shell crossing occurs. 

When following the respective solutions into the non¬ 
linear regime we face the problem that the gauge trans¬ 
formation becomes highly nontrivial. This makes it diffi¬ 
cult to compare quantities like the density or metric com¬ 
ponents directly, because they are gauge-dependent. A 
good way to proceed is to compute observables, which are 
by definition gauge-independent. In the following section 
we will discuss some examples in detail. 


C. Observables 

1. Redshift of radially in-falling source 

The first observable we will study is the redshift of 
a source of light that is moving with the flow of par¬ 
ticles surrounding it. We place this source of light at 
an initial radius r\ from the center and an observer at 
r-i, the boundary where the inhomogeneous LTB patch 
is matched to FLRW. For this example, the source and 
observer are along the same radial line. The source con¬ 
stantly emits photons at a fixed energy given in the rest 
frame of the source. As the simulation progresses, we 
propagate these photons through the simulation volume 
until they reach the observer. There they are detected 
and we calculate their observed redshift, which is defined 
as: 


to the momentum p of the source (observer) particle in 
the limit of weak fields as: 


(1-I-'k)— ra^ + p^ — \ 


(17) 


To propagate a photon through the simulation, we use 
the null condition (ds^ = 0). We can actually find a 
fixed relation between dr/dr and dp/dr, a consequence 
of the fact that a photon always travels at the speed of 
light: 


1 -I- 2^' -I- 2$ = 



which gives: 


dr 

dr 


±(l + «'-k$) 






(18) 


(19) 


On every step, the photon’s energy can be evaluated by 
integrating the time component of the geodesic equation: 


dfcO 

dr 




dr 

dr 


+ 2n 


k° = 0 


( 20 ) 


The results for this observable are shown in figure 4. 
As can be seen, the simulation agrees well with the LTB 
predictions. In fact, for this plot the leading source of de¬ 
viation actually comes from imprecise matching of initial 
conditions, which could only be improved by performing 
that matching at a higher order of perturbation theory. 


2. Lensing of non-radial rays 


Another observable we can analyse is the deflection of a 
ray that propagates not only in the radial, but also in an 
angular direction. The trajectory of such a ray is lensed 
by the gravitational potentials. Spherical symmetry en¬ 
sures that the trajectory of a light ray will be planar, so 
we only have to consider the radial direction and one an¬ 
gular direction. By setting i? = 7r/2 or = 0, one can 
derive the two equations that determine that path of the 
photon from the geodesic equation: 


d^r 

dr^ 


2(4-,^ -k4>,r) 



(4',^-k4>,r ) 


dr 

dr 


dp 

dr 


r -|- (<i>,r -|-'I',j.) — 0 


(21a) 


d^p 

dr^ 


-k(4>,r-4',r 





dr dp ^ 
dr dr 


(21b) 


1 T -^obs 


)|src 

{gf_„ykf^u‘')\ohs ’ 


(16) 


where the product of the photon’s 4-momentum and 
the 4-velocity of the source (observer) can be related 


Varying the initial angle at which photons enter the per¬ 
turbed region and observing the deflection angle (the an¬ 
gle by which its outgoing trajectory differs from the in¬ 
coming one), we get a gauge independent probe of the 
underlying potentials. 
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Figure 2: (Color online) Top left: the evolution of the density profile of a spherically symmetric compensated tophat pertur¬ 
bation. Different-coloured lines correspond to outputs at different times in the simulation, parametrised by the background 
redshift «flrw- The last output at ^flrw = 0.5 happens just before the collapse occurs. Density plots exhibit some discreteness 
noise, which is caused by having a finite number of particles. Bottom left: the momentum of the shells moving inwards. Top 
right: evolution of the underlying scalar metric perturbation $. This profile is continuous even where the density has a step. 
Bottom right: Evolution of the difference of the two potentials x = ~ This is a purely relativistic quantity and does not 

exist in a Newtonian setup. The magnitude of x is ~ ■ Parameters of the simulation were: size of the box: 20 Mpc//i, initial 

radii of top-hat overdensity and compensated region, respectively: ri = 6 Mpc/h, r 2 = 18 Mpc//i, initial density contrast of 
the overdensity: 5 = 1/200, initial redshift: Zin = 500. 


In Figure 5 we show the trajectories of photons that 
propagate through the simulation volume. We tracked 
200 photons, entering at different angles a. The photons 
that experience the most lensing are those that pass near 
the edge of the overdensity. A photon that enters at 
a = 0 and passes through the centre of the overdensity 
is not lensed, since its path is radial. Likewise, a photon 
that enters at a = 7^12 spends too little time inside the 
non homogeneous region to change its path substantially. 
In Figure 6 we show the deflection angle as a function of 
a. 


III. ANGULAR MOMENTUM 


Using the formalism presented so far, any initial over¬ 
density would collapse to a singularity within a finite 
amount of time. For realistic cosmological structures this 
does not happen due to the process of virialisation during 
which the initial potential energy of a large scale density 
fluctuation is partially converted into radial and angu¬ 
lar kinetic energies of individual particles. Most impor¬ 
tantly, the angular momentum thus generated in individ¬ 
ual particles causes these particles to miss the centre of a 
collapsing structure. This then avoids the production of 
densities large enough to cause singularities to arise. Es- 
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Figure 3: (Color online) The same set of plots as for figure 2, this time for an evolving underdensity. The last two outputs 
exhibit a shell-crossing, which can be seen in the density and momentum portraits. This non-linear feature can only be modelled 
with an N-body simulation. The parameters used here were: size of the box: 120 Mpc//i, ri = 36 Mpc/h, r 2 = 84 Mpc//i, 
S = -1/40, Zin = 500. 


sentially, the produced angular momentum provides an 
effective pressure term that resists the collapse. 

We show in this section how we can model this pres¬ 
sure by adding angular momentum to the particles in 
our simulation box, without losing spherical symmetry. 
The downside to the method we present is that there can 
be no exchange of angular momentum between particles. 
This is because adding such an effect would necessarily 
require some degree of deviation from spherical symme¬ 
try. Unfortunately, this limits how far we can model the 
true virialisation process; nevertheless, as we show, we 
can still set initial conditions that produce stable, bound, 
spherical structures. 

We can now imagine the shells to be made up of in¬ 
finitesimal point-like particles, evenly distributed over 
the sphere. Apart from the radial momentum, which 
is the same for all infinitesimal point-like particles on a 
given sphere, we can assign each of those particles an an¬ 


gular momentum in a particular direction, in such a man¬ 
ner that once we perform the average over the momenta 
of all infinitesimal particles on a given sphere, there is 
no preferred direction of angular momentum. One can 
imagine that for every infinitesimal particle with some 
angular momentum, there is another particle on a tra¬ 
jectory in the same plane, but traveling in the opposite 
direction. Although we can not observe any non-radial 
motion of spherical shells, their radial motion is neverthe¬ 
less affected. This is because the equation that governs 
the propagation of particles in the radial direction in¬ 
volves a pressure-like term that depends on the angular 
momentum of the particle. 

We start again with the geodesic equation 


, dx^n) 

ds 


ds^ 


ds 


= 0 


( 22 ) 


Calculating all the Christoffel symbols we have three 




































Figure 4: (Color online) Top: the redshift of light, emitted 
by an in-falling source particle, as a function of the back¬ 
ground redshift zflrw at the time when the light is detected 
by the observer. Initially, the observed redshift is decreas¬ 
ing, which is something we would expect for two particles co¬ 
moving with the background in a matter dominated universe. 
Later, as the collapsing structure evolves, the velocity of the 
in-falling source becomes the dominant contribution and the 
redshift starts increasing again. Bottom: relative error be¬ 
tween our relativistic simulation and the LTB solution. The 
error is mainly due to the first-order matching of the initial 
conditions. The error increases as the collapse evolves. This 
is because the collapse time itself receives a first-order correc¬ 
tion. Therefore, the divergence in observed redshift happens 
at slightly offset times, resulting in the error blowing up. For 
this plot, the same parameters as the ones in Figure 2 were 
used. 



Figure 5: Schematic representation of the trajectories of 
lensed photons in an LTB geometry. The rays enter the LTB 
region at varying angles on the right end of the plot. Along 
their way through the simulation volume their trajectories are 
deflected due to the underlying overdensity. 



equations for evolution of coordinates: 








- 



dr dr 
ds ds 
\ 2 



(23a) 


Figure 6: (Color online) Tensing of photons by an underlying 
overdensity. We plot the LTB predictions (solid lines) and the 
results from our simulation (dashed lines on top of them) at 
four selected redshifts. The parameters used here are the 
same as the ones in Figure 2. 


dV 

ds^ 


-2{n 




dTd^_2 
ds ds 



dr dip 
ds ds 


(23c) 


d^r 

ds^ 


= -T'. 





- 2 {n-d>^r) 

-I- r (1 — r<i)_r) 


dr dr 
ds ds 



Here we only kept one angular coordinate. Since the 
potentials are spherically symmetric, each infinitesimal 
point-like particle will move in a two-dimensional plane 
and so the reference frame can always be rotated so that 
the direction of the angular momentum is perpendicular 
to this plane. We see that the last equation (23c) can be 
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Figure 7: (Color online) Space-time diagrams showing the radial trajectories of shells initially arranged as a uniform ball with 
zero radial motion. Different colors label different shells, black labelling the outermost one. From left to right, the angular 
momentum of the shell particles was set to L = 0 where we get collapse; L = O.SLk where a violent re-configuration to a 
new, much more compact state happens; L = Lk at which point the behaviour is almost Keplerian, with tiny perturbations 
caused by relativistic corrections - stability is guaranteed for a very long time scale; L = 1.1 Lk with quasi-stable oscillations 
with relatively long timescale for chaotic behavior; and L = 1.4 Lk which corresponds to an unbound state. Flere, Lk is the 
(r-dependent) value of L which corresponds to a circular Keplerian orbit in Newtonian theory. The radius is plotted in units 
of rs, the Schwarzschild radius of the ball; initially the ball has a radius of 50rs. The time coordinate is plotted in units of tc, 
the collapse time of the irrotational ball. 


integrated analytically with respect to ds once, to give 

^ = -^(1 + 2 $) ( 24 ) 

ds (ar)^ 

where L is the constant of integration. We can express 
the evolution of coordinates r and ip with respect to co¬ 
ordinate time dr instead of eigentime ds using a simple 


trick: 

d^r _ /d^r 

dr^ yds^ 

And equivalently for ip. With this we can combine equa¬ 
tions (23b) and (23c) with (23a) to express: 


d^T dr 
ds^ dr 


dr \ 
ds / 


-2 


(25) 


d^r 

dr^ 


dr 

{-U + 2^ r + t) - 

dr 


(4>,r -b 2'i>,r) 





-b (-^(1 - 2 $ - 24 -) - 4 >,^) 


2 


dr 

dr 


(26a) 
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2 rx T \ 

- h 2$ y. + 2^ j. I ——-— 

r / dr dr 


(26b) 


r 


In addition we have the “mass-shell condition” for par¬ 
ticles with mass: 


dx^^ dx'' 

which is expressed, using our metric, as:^ 


(27) 


-fo2r2(^) (1-2$) = -! (28) 

The momentum of the particles now has a radial compo¬ 
nent. 


Pr = 


m(l-d>)f 


+ 2VI/ - (1 - 2$) - (1 - 2ci>)r2 


and an angular one, 
Pv = 


(29) 


m(l - ^')r^ 


+ 2VI/ - (1 - 2$) (^)" - (1 - 2<i>)r2 

(30) 

Using the mass-shell condition (28), the definition of 
the conserved quantity L (24), and (23a), we can extract 
the angular velocity of infinitesimal particles. 


dp L 
dr ar^ ^ 


1-h 4$-1-24-- (1-f 2$) (^) 
1+Wt) 


(31) 


and eliminate it from the momenta equations: 


Pr = m 


dr(l-4>) / a?r^ + L'^{l + 2^) 


d'T ar Y 1 + 24--(1-2$) ()^)' 


and 


Pip = 


mL(l -I- $) 


ar 


(32) 


(33) 


Note that it is L which is the conserved quantity, not p^p, 
within our framework. With this setup, we can finally 


® Again, we have set = 7r/2, so di9/ds = 0. 


express the equation for evolution of shell particles (23b) 
as a system of two first-order equations: 


^ = (^,r - 'H)Pr - ^m? + pI + p%dl^r 

pl{l + ^ + d!) 


- $ 


and 


dr 

dr 


Pr 


+pI+ Pp 


+Pr+ pI 


;(l-h$-h$). 


(34) 


(35) 


These two equations describe the movement of shell par¬ 
ticles and need to be solved numerically. A sanity check 
verifies that if we set to zero, we recover the mo¬ 
mentum evolution equation in the case without angular 
momentum, (12). 

With our new definitions of Pr and p^, the energy- 
momentum tensor expresses as: 

.0 1 + 3$ V- 


^0 = {\/^"+Pr+Pl) (36) 


and 



TT _ 1 + 3$ 
“ 47rr2a3 


The angular motion also gives rise to a transverse 
Doppler effect. This can be seen, e.g., from eq. (17) being 
modified as 


Qfi.uk^u'' = -k°a{l -I- $) — [a/t 

m L V 


■Pr 


■P% 


-Pr 


(38) 

In figure 7 we have plotted the radial trajectories of 
shells within balls that have uniform density, but non¬ 
zero and non-uniform angular momentum. Each panel 
corresponds to the same initial density state, but differ¬ 
ent initial states of angular momentum. As can be seen, 
our code is able to describe balls that collapse to a point; 
balls that first begin to collapse under gravity but then 
stabilise; balls where the effective pressure due to an¬ 
gular momentum perfectly balances gravitational attrac¬ 
tion; balls that first expand due to effective pressure, but 
then stabilise; and finally, balls which are blown apart by 
pressure. Only in the first and last situations respectively 
will our code certainly break down, and even then it will 
survive until the weak field limit breaks down, or a par¬ 
ticle leaves the box. In these plots, the balls begin only 
50 times larger than their Schwarzschild radii, therefore 
relativistic effects will not be negligible. 
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IV. CONCLUSION 

We have presented an N-body framework for spheri¬ 
cally symmetric solutions valid in the weak-field regime of 
general relativity. We have primarily applied this frame¬ 
work in a cosmological context, expanding around an 
FLRW metric; however nothing forbids the application 
of the framework to other contexts. Spherical symmetry 
was imposed in order to obtain an economic setup for nu¬ 
merical studies. We compared our code against two types 
of known exact solutions, Schwarzschild and Lemaitre- 
Tolman-Bondi, and found good agreement. However, our 
scheme is suitable also for setups where no exact solution 
is known, for instance when the fluid description of mat¬ 
ter is not valid. 

Furthermore, we demonstrated that the relativistic po¬ 
tentials are computed more accurately than in a Newton¬ 
ian scheme. This feature will be useful for the study of 
models which have exotic sources of stress-energy per¬ 
turbations, such as dynamical dark energy or modified 
gravity. On a related note, we stress that our scheme 
does not make any assumptions about the nature of per¬ 
turbations apart from the requirement that they give rise 
to weak gravitational fields only. This assumption breaks 
down, e.g., if a black hole is formed. 

In order to avoid the collapse of an overdensity into 
a black hole, we have introduced a method to create a 
stable bound structure supported by angular momentum. 
Such configurations may, in some sense, be more realistic 
proxies for cosmic structures such as galaxy clusters, and 
can therefore be useful laboratories for studying gravity 
at these scales. They may also be useful for the study of 
weakly relativistic, compact bound objects that can form 
in the early universe, such as ultra compact mini-haloes, 
or for the early stages of the formation of primordial black 
holes. 
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Appendix A: Linear Relation between LTB Solution 
and Longitudinal Gauge 

In this appendix we give the linear gauge transforma¬ 
tions which we use to set up initial conditions correspond¬ 
ing to a given LTB solution. These relations are valid to 


the extent that the matching is done at a time when lin¬ 
ear perturbation theory can be applied. LTB solutions 
which do not allow a perturbative description at any time 
are more difficult to translate into our gauge, and there 
may be cases where it is impossible. We will not consider 
solutions of this type in this paper. 

Our starting point is the LTB line element given in 
eq. (15). At the initial time tin where we will do the 
matching, we will rescale the coordinate r such that, at 
that time, i?(tin,r) = a(tin)r. The (time independent) 
gravitational mass function can then simply be obtained 
as 

2M(r) = 87 rGa^(tin) f f^p(tin,f)dr 

Jo 

= [ f2(l + ,5(tin,f))dr. (Al) 

Jo 

In this work we are interested in setups where the metric 
is FLRW everywhere except for a finite spherical region. 
This can be achieved by choosing a “compensated” den¬ 
sity profile such that the mass function becomes the one 
of FLRW, M(r) = /2, at the boundary of the re¬ 

gion^. We will choose a particularly simple profile, given 

by 

{ 5i r < ri 

82 ri < r < r 2 , (A2) 

0 r > r 2 

where 82 = —(5irf/ (rf — rj) gives the correct matching 
to FLRW as can be seen by inspecting the resulting mass 
function. 

Next we will use the parametric expression for the ex¬ 
act LTB solution in order to determine the metric func¬ 
tion E{r). Let us consider E{r) < 0, the opposite case is 
analogous. The parametric expressions are 

R{t,r) = -^^||^(l-cos? 7 ) , (A3a) 

(^Sb) 

While the parameter 77 cannot be eliminated in closed 
form, it is possible to do so perturbatively for small 77 , 
i.e. at early time. We And that 

E{r) = - ^(tin)a^ (tin)/(?-), (A4) 

where 

(81 if r < ri 

f{i")=l 82 +^{ 8 i- 82 ) if ri < r < r 2 • (AS) 

I 0 if r > r 2 


Note that is independent of time in a matter dominated 

universe. 





12 


In the linear regime we also have which implies 

(A6) 

Let us now write the line element in a convenient pertur¬ 
bative notation, 


(A8.) 

y{t,r) = b{t,r) + rb^r{t,r) — E{r). (A8b) 


r) = a{t)r 


a{tl) 


R\t,r) 

[R^r{t,r)f 

1 -h 2E{r) 


a^{t)r'^ [1 -h 2 b{t,r)] , 
a^{t) [l + 2y{t,r)\ , 


(A7a) We can now work out the linear gauge transformation 
from synchronous to longitudinal gauge (see also [28]). 
(A7b) A straightforward calculation shows 


dr 

d 


r 

r 


-b^r-r{t, r) + H{t)a^{t) [y,t(t, r) - r)] - ^ [y{t, r) - b{t, r)] -k r), (A9a) 

2H{t)a^{t) [6,t(t, r) - y^t{t, r)] + a^{t) [6,tt(t, r) - y._tt{t, r)] . (A9b) 


Since the combination is independent of time in a 

matter dominated universe we can evaluate it at t = tin 
and find 


dr 


-'^,r{t,r) 

r 


= r 


dr 


-^,rit,r) 


= -H\u^)a\U^)rf{r), (AlO) 


which can be integrated twice to obtain and dt. The 
constants of integration should be chosen such that we 
can match smoothly to FLRW at r = r 2 . In other words, 
we require ^\r=r2 = '^\r=r2 = ^,r\r=r2 = ^ ,r\r=r2 = 0 - 
The corresponding solutions are 

$(t,r) = ^'(t,r) = ii7^(ti„)a^(ti„) / r/(r)df. (All) 


Appendix B: Initial Particle Data 

Given a linear solution for <&, d' which specifies the 
initial conditions, we can work out the initial particle 
configuration. To this end, we linearize eq. (2), 


To see this, take the continuum limit of the particle sum. 


^0 _ 

-^0 — 


l-k3$ 

47rr^a^ 


^TO(„)5(r-r(„)) 


[ fK))Sir - r(„))dr°„) , (B3) 


with a distribution function f{r°^) oc (r°„^)^ correspond¬ 
ing to the homogeneous distribution. Next, change the 
integration variable to r(„) to obtain eq. (B2). Insert¬ 
ing into eq. (Bl) and using eq. (3), the solution for the 
displacement is found to be 

Here the constant of integration is fixed by requiring reg¬ 
ularity at the origin, which implies 5r|r=o = 0. 

The initial particle velocities can be obtained simply 
by taking the time derivative of above equation. 


dr 

dr 


tin 


Sr. 



(B5) 




^ - 377 ^$ = -AirGa^STS , 

r ' 


(Bl) 


where we used R' = —77^/2 in a matter dominated uni¬ 
verse. 


where we have used that <i> t = y; = 0 at linear order. 

The aim of this section is to construct a linear dis¬ 
placement field 5r(r) which specifies the initial particle 
positions r(„)(tin) = -|- 5r(r°^j) as infinitesimal dis¬ 
placements from a homogeneous distribution Ex¬ 

panding eq. (8) to linear order we find 

ST^ = ^3$ - Sr^r - j . (B2) 


Appendix C: Particle-Mesh Scheme for Spherical 
Coordinates 

Standard particle-mesh schemes [29] usually employ a 
Cartesian mesh which means that a few modifications are 
required in order to make them fit for our purpose. Our 
mesh will have a uniform resolution in r, the radial co¬ 
ordinate, meaning that the volume of the cells increases 
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as one moves outwards from the center. The mass reso¬ 
lution can be set independently by changing the number 
of particles per cell - a number which can also depend on 
radius and should be chosen according to the problem at 
hand. 

The stress-energy tensor on the grid is computed by 
means of a particle-to-mesh projection. It is constructed 
by smearing out each particle over a finite radial interval 
and then determining the fraction of its mass within each 
cell. Explicitly, we replace 


^ w{r, - r(„)) 

47 ^ 7-2 y. 


(Cl) 


where denotes the center of the ith cell, Vi is the cell’s 
volume, r(n) is the position of the nth particle, and w is 
a weight function which depends on the smearing. We 
use triangular-shaped particles where 


w{ri 



ri)A(r-r(„))dr, 


(C2) 


with 


n(Ar) 


1 if - ^ < Ar < f 
0 otherwise 


(C3) 


and 


for the cell at the origin. Contributions which would be 
projected at negative radius are simply folded back onto 
the positive axis. 

Next, we want to construct a particle distribution 
which would correspond to a homogeneous universe. A 
perturbed distribution can then be obtained by acting 
with an infinitesimal displacement as explained in ap¬ 
pendix B. Our homogeneous distribution will be con¬ 
structed in a way as to minimize discretization issues. 
For simplicity, let us discuss a setup where we have one 
particle per grid cell (mass resolution can be increased 
by subdividing particles). If we choose initial particle 
positions as r(„) = {n+ |)dr, with n = 0,1,..., one can 
recursively construct the mass assignment for each par¬ 
ticle which would lead to an exactly uniform projected 
density under the above projection method: 

m(„) = 47r + n + dr^a^p (C7) 

Here, p = —Tq is the homogeneous density of the back¬ 
ground FLRW model. Evidently, at very large radius 
r(„) dr this expression asymptotes to the correct con¬ 
tinuum limit TO(„) = 47rr^^jdra^p. The corrections which 
come in at small radii are chosen to compensate for dis¬ 
cretization effects, such that the projected density re¬ 
mains exactly homogeneous. 


A(Ar) 


1 -k ^ if - dr < Ar < 0 
1 -^ ifO<Ar<dr 
0 otherwise 


(C4) 



w{ri 


- Hn)) 

V 


^(n) 


P 


Vz. 


(C8) 


Here and in the following, dr denotes the grid unit. Picto- 
rially, □ characterizes the footprint of the cell (an interval 
of width dr centered at r = r^), whereas A characterizes 
the shape of the particles (its mass distribution along the 
radial coordinate). 

Our grid cells are centered at = fdr, with i = 
0,1,.... The volume of the Ah cell is computed as 




47r 

T 





3 


= 47r 




(C5) 


for i > 0, and 




(C6) 


The weight function w can also be used in order to in¬ 
terpolate grid-based quantities (fields, gradients of fields 
etc.) to the positions of the particles. This is neces¬ 
sary for the integration of the geodesic equations. For 
instance, in order to interpolate 4>, we would define 


Note that the sum is now taken over the grid points. 
Similarly, we can interpolate a gradient as 


= X — 2 ^—-w^(a + -dr - r(„)) 


(CIO) 

based on a standard one-sided two-point gradient which 
naturally sits at half-integer grid units. 
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